      SUBROUTINE CALSTEPD  
C  
C CHANGE RECORD  
C **  SUBROUTINE CALSTEP ESTIMATE THE CURRENT MAXIMUM TIME STEP SIZE  
C **  FORM LINEAR STABILITY CRITERIA AND A FACTOR OF SAFETY  
C  
      USE GLOBAL  

      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::DTL1  
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::DTL2  
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::DTL3  
      REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)::QSUBINN  
      REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)::QSUBOUT  
      IF(.NOT.ALLOCATED(DTL1))THEN
        ALLOCATE(DTL1(LCM))
        ALLOCATE(DTL2(LCM))  
        ALLOCATE(DTL3(LCM))  
        ALLOCATE(QSUBINN(LCM,KCM))  
        ALLOCATE(QSUBOUT(LCM,KCM))  
        DTL1=0.0 
        DTL2=0.0 
        DTL3=0.0 
        QSUBINN=0.0 
        QSUBOUT=0.0 
      ENDIF        
C
C**********************************************************************C
C
      ITRNTMP=0
      DO NX=1,7
        ITRNTMP=ITRNTMP+ISTRAN(NX)
      ENDDO
C
      IF(N.LE.0)DTDYN=0.0
C
      DTMIN=DT
      DTMAX=TIDALP
C
      DO L=2,LA
        DTL1(L)=DTMAX
        DTL2(L)=DTMAX
        DTL3(L)=DTMAX
      ENDDO
C
C **  DETERMINE SOURCE/SINKS FOR SUBGRID SCALE CHANNEL EXCHANGES
C
      DO K=1,KC
        DO L=2,LA
          QSUBOUT(L,K)=0.0
          QSUBINN(L,K)=0.0
        ENDDO
      ENDDO
C
      IF(MDCHH.GE.1)THEN
        DO K=1,KC
        DO NMD=1,MDCHH
          LMDCHHT=LMDCHH(NMD)
          LMDCHUT=LMDCHU(NMD)
          LMDCHVT=LMDCHV(NMD)
          IF(MDCHTYP(NMD).EQ.1)THEN
            QUKTMP=QCHANU(NMD)*DZC(K)
            QVKTMP=0.
          ENDIF
          IF(MDCHTYP(NMD).EQ.2)THEN
            QVKTMP=QCHANV(NMD)*DZC(K)
            QUKTMP=0.
          ENDIF
          IF(MDCHTYP(NMD).EQ.3)THEN
            QUKTMP=QCHANU(NMD)*DZC(K)
            QVKTMP=QCHANV(NMD)*DZC(K)
          ENDIF
          QSUBOUT(LMDCHHT,K)=QSUBOUT(LMDCHHT,K)
     &               +MIN(QUKTMP,0.)
     &               +MIN(QVKTMP,0.)
          QSUBINN(LMDCHHT,K)=QSUBINN(LMDCHHT,K)
     &               +MAX(QUKTMP,0.)
     &               +MAX(QVKTMP,0.)
          QSUBOUT(LMDCHUT,K)=QSUBOUT(LMDCHUT,K)
     &               -MAX(QUKTMP,0.)
          QSUBINN(LMDCHUT,K)=QSUBINN(LMDCHUT,K)
     &               -MIN(QUKTMP,0.)
          QSUBOUT(LMDCHVT,K)=QSUBOUT(LMDCHVT,K)
     &               -MAX(QVKTMP,0.)
          QSUBINN(LMDCHVT,K)=QSUBINN(LMDCHVT,K)
     &               -MIN(QVKTMP,0.)
        ENDDO
        ENDDO
      ENDIF
C
C **  METHOD 1: UPWIND DIFF IN MOMENTUM EQUATIONS
C
      DO K=1,KC
        DO L=2,LA
          IF(LMASKDRY(L))THEN
          LE=LEAST(L)
          LN=LNC(L)
          LS=LSC(L)
          KM=K-1
          VATUUU=0.25*(V(L,K)+V(LWEST(L),K)+V(LN,K)+V(LN-1,K))
          TMPUUU=ABS(U(L,K)/DXU(L))+ABS(VATUUU/DYU(L))
          DTTMP=1./(TMPUUU+1.0E-18)
          DTL1(L)=MIN(DTL1(L),DTTMP)
          UATVVV=0.25*(U(L,K)+U(LS,K)+V(LEAST(L),K)+V(LS+1,K))
          TMPVVV=ABS(V(L,K)/DYV(L))+ABS(UATVVV/DXV(L))
          DTTMP=1./(TMPVVV+1.0E-18)
          DTL1(L)=MIN(DTL1(L),DTTMP)
          UEAST=ABS(U(L,K))
          UWEST=ABS(U(LEAST(L),K))
          VSOUTH=ABS(V(L,K))
          VNORTH=ABS(V(LN,K))
          TMPVVV=MAX(VSOUTH,VNORTH)
          TMPUUU=MAX(UEAST,UWEST)
          TMPVAL=TMPUUU/DXP(L)+TMPVVV/DYP(L)
          DTTMP=1./(TMPVAL+1.0E-18)
          DTL1(L)=MIN(DTL1(L),DTTMP)
          ENDIF
        ENDDO
      ENDDO
C
C **  METHOD 2: POSITIVITY OF ADVECTED MATERIAL, DTL2
C
      IF(ITRNTMP.GE.1)THEN
      DO K=1,KC
        DO L=2,LA
          !IF(LMASKDRY(L))THEN
          ! *** DSLLC SINGLE LINE (TO ADDRESS BED CHANGE, IMORPH UPDATES H1P TO HP
          IF(HP(L).GT.HDRY.AND.H1P(L).GT.HDRY)THEN  
          LE=LEAST(L)
          LN=LNC(L)
          LS=LSC(L)
          KM=K-1
          TOP=DZC(K)*H1P(L)*DXYP(L)
          QXPLUS=UHDY2(LE,K)*DZC(K)
          QXPLUS=MAX(QXPLUS,0.0)
          QYPLUS=VHDX2(LN,K)*DZC(K)
          QYPLUS=MAX(QYPLUS,0.0)
          QZPLUS=W2(L,K)*DXYP(L)
          QZPLUS=MAX(QZPLUS,0.0)
          QXMINS=UHDY2(L,K)*DZC(K)
          QXMINS=-MIN(QXMINS,0.0)
          QYMINS=VHDX2(L,K)*DZC(K)
          QYMINS=-MIN(QYMINS,0.0)
          QZMINS=W2(L,KM)*DXYP(L)
          QZMINS=-MIN(QZMINS,0.0)
          QTOTAL=QSUM(L,K)+QSUBOUT(L,K)+QSUBINN(L,K)
          QSRC=-MIN(QTOTAL,0.0)
          BOT=QXPLUS+QYPLUS+QZPLUS+QXMINS+QYMINS+QZMINS+QSRC
          IF(BOT.GT.0.0)THEN
            DTTMP=TOP/BOT
            DTL2(L)=MIN(DTL2(L),DTTMP)
            IF(DTTMP.LT.0.0)THEN
              WRITE(6,880)IL(L),JL(L),K,TOP,QXPLUS,QYPLUS,QZPLUS,
     &                    QXMINS,QYMINS,QZMINS,QSRC
              WRITE(8,880)IL(L),JL(L),K,TOP,QXPLUS,QYPLUS,QZPLUS,
     &                    QXMINS,QYMINS,QZMINS,QSRC
            ENDIF
          ENDIF
          ENDIF
        ENDDO
      ENDDO
      ENDIF
C
C **  METHOD 3: implicit BOTTOM FRICTION AND ROTATIONAL ACCELERATION DAMPING
C
      DO L=2,LA
        IF(LMASKDRY(L))THEN
          TMPVAL=SUB(L)+SUB(LEAST(L))+SVB(L)+SVB(LNC(L))
          IF(TMPVAL.LT.0.5)THEN
            LN=LNC(L)
            TAUBC=QQ(L,0)/CTURB2
            UCTR=0.5*(U(L,1)+U(LEAST(L),1))
            VCTR=0.5*(V(L,1)+V(LN,1))
            UHMAG=HP(L)*SQRT(UCTR*UCTR+VCTR*VCTR)
            IF(UHMAG.GT.0.0)THEN
              FRIFRE=TAUBC/UHMAG
              FRIFRE2=FRIFRE*FRIFRE
              ACACTMP=(CAC(L,KC)*HPI(L)*DXYIP(L))**2
              IF(ACACTMP.GT.FRIFRE2)THEN
                DTTMP=2.*FRIFRE/(ACACTMP-FRIFRE2)
                DTL3(L)=MIN(DTL3(L),DTTMP)
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDDO
C
C **  CHOOSE THE MINIMUM OF THE THREE METHODS
C
      DTL1MN=2.*DTMAX
      DTL2MN=2.*DTMAX
      DTL3MN=2.*DTMAX
      DTTMP=2.*DTMAX
      DO L=2,LA
        IF(LMASKDRY(L))THEN
          IF(DTL1MN.GT.DTL1(L))THEN
            DTL1MN=DTL1(L)
            L1LOC=L
          ENDIF
          IF(DTL2MN.GT.DTL2(L))THEN
            DTL2MN=DTL2(L)
            L2LOC=L
          ENDIF
          IF(DTL3MN.GT.DTL3(L))THEN
            DTL3MN=DTL3(L)
            L3LOC=L
          ENDIF
        ENDIF
      ENDDO
C
C *** DSLLC BEGIN BLOCK
C
C **  FIND MINIMUM & APPLY A SAFETY FACTOR
C
      DTL1MN=DTL1MN*DTSSFAC
      IF(DTTMP.GT.DTL1MN)THEN
        DTTMP=DTL1MN
        DTCOMP=DTTMP/DTSSFAC
        LLOC=L1LOC
      ENDIF
      DTL2MN=DTL2MN*DTSSFAC
      IF(DTTMP.GT.DTL2MN)THEN
        DTTMP=DTL2MN
        DTCOMP=DTTMP/DTSSFAC
        LLOC=L2LOC
      ENDIF
      DTL3MN=DTL3MN*DTSSFAC
      IF(DTTMP.GT.DTL3MN)THEN
        DTTMP=DTL3MN
        DTCOMP=DTTMP/DTSSFAC
        LLOC=L3LOC
      ENDIF
      LMINSTEP=LLOC
C *** DSLLC END BLOCK
C
C **  CHECK IF CURVATURE INSTABILITY IS CONTROLLED
C
c      CACDTMX=-1000.
c      DO L=2,LA
c        IF(LMASKDRY(L))THEN
c          CACTMP=ABS(DTTMP*CAC(L,KC)*HPI(L)*DXYIP(L))
c          CACDTMX=MAX(CACDTMX,CACTMP)
c        ENDIF
c      ENDDO
c      CACAMP=SQRT(1.+CACDTMX*CACDTMX)
C
      TIMEDAY=TIMESEC/86400.
      IF(DTCOMP.LT.DTMIN)THEN   ! *** DSLLC SINGLE LINE
        WRITE(8,800)TIMEDAY,DTTMP,DTMIN,IL(LLOC),JL(LLOC)
        WRITE(6,800)TIMEDAY,DTTMP,DTMIN,IL(LLOC),JL(LLOC)
        WRITE(8,801)IL(L1LOC),JL(L1LOC),DTL1MN
        WRITE(6,801)IL(L1LOC),JL(L1LOC),DTL1MN
        WRITE(8,802)IL(L2LOC),JL(L2LOC),DTL2MN
        WRITE(6,802)IL(L2LOC),JL(L2LOC),DTL2MN
        WRITE(8,803)IL(L3LOC),JL(L3LOC),DTL3MN
        WRITE(6,803)IL(L3LOC),JL(L3LOC),DTL3MN
        DTTMP=DTMIN
C *** DSLLC BEGIN BLOCK
      ELSEIF(DTTMP.LT.DTMIN)THEN
        DTWARN=DTTMP
        DTTMP=DTMIN
C *** DSLLC END BLOCK
      ELSE
        TMPVAL=DTTMP/DTMIN
        ITMPR=NINT(TMPVAL)
        RTMPR=FLOAT(ITMPR)
        IF(RTMPR.LT.TMPVAL)THEN
          DTTMP=RTMPR*DTMIN
        ELSE
          DTTMP=(RTMPR-1.)*DTMIN
        ENDIF
      ENDIF
C
C **  SET TO MINIMUM TIME STEP ON STARTUP
C
      IF(N.EQ.0)DTTMP=DTMIN
C
C **  RESTRICT INCREASE IN TIME STEP TO DTMIN
C
C      DTDYN2=2.*DTDYN
C      IF(DTTMP.GT.DTDYN2)THEN
C        DTTMP=DTDYN2
C      ENDIF
C
      DTDYNP=DTDYN+DTMIN
      IF(DTTMP.GT.DTDYNP)THEN
        DTTMP=DTDYNP
      ENDIF
C
      DTDYN=DTTMP
C
C **  SET INCREMENTAL INCREASE IN OUTPUT COUNTER
C
      NINCRMT=NINT(DTDYN/DTMIN)
C
C **  ADJUST INCREMENT FOR N TO LAND EVENLY ON NTSPTC
C
      RTCTMP=FLOAT(N)/FLOAT(NTSPTC)
      NTCTMP=RTCTMP
      NTMP=(1+NTCTMP)*NTSPTC-N
      IF(NINCRMT.GT.NTMP)THEN
        NINCRMT=NTMP
        DTDYN=FLOAT(NTMP)*DTMIN
      ENDIF
C
  100 FORMAT(5I5,5F12.5,E13.5)     
  101 FORMAT(3I5,E13.5)     
  800 FORMAT('  TIME,DTDYN,DTMIN,I,J = ',F12.5,2E12.4,2I7)     
  801 FORMAT('  MOM  ADV,I,J,DTM = ',2I5,E13.4)     
  802 FORMAT('  MASS ADV,I,J,DTM = ',2I5,E13.4)     
  803 FORMAT('  CURV ACC,I,J,DTM = ',2I5,E13.4) 
  880 FORMAT(3I5,8E13.4) 
 8899 FORMAT(' DT3 ERROR ',2I5,6E13.5)   
C
C**********************************************************************C
C
      RETURN
      END
